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Introduction 

Approximately  2.5  million  Americans  live  with  epilepsy  and  epilepsy-related  deficits  today,  more 
than  disabled  by  Parkinson  disease  or  brain  tumors.  The  impact  of  epilepsy  in  the  US  is 
significant  with  a  total  cost  to  the  nation  for  seizures  and  epilepsy  of  approximately  $12.5  billion. 
Epilepsy  consists  of  more  than  40  clinical  syndromes  affecting  40  million  people  worldwide. 
Approximately  25  percent  of  individuals  receiving  antiepileptic  medication  have  inadequate 
seizure  control;  however,  80%  individuals  with  medication  resistant  epilepsy  might  be  cured 
through  surgery  if  one  were  able  to  precisely  localize  the  seizure  focus.  The  proposed  research 
will  significantly  advance  our  ability  to  localize  such  foci,  and  thereby  offer  curative  epilepsy 
surgery  for  this  devastating  disease.  Photoacoustic  tomography  (PAT)  uniquely  combines  the 
high  contrast  advantage  of  optical  imaging  and  the  high  resolution  advantage  of  ultrasound 
imaging  in  a  single  modality.  In  addition  to  high  resolution  structural  information,  the  proposed 
PAT  is  also  able  to  provide  functional  information  that  are  strongly  correlated  with  regional  or 
focal  seizure  activity,  including  blood  volume  and  blood  oxygenation  because  of  the  high 
sensitivity  of  optical  contrast  to  oxyhemoglobin  and  deoxyhemoglobin  concentrations.  The 
hypothesis  of  the  proposed  research  is  that  PAT  offers  the  possibility  to  non-invasively  track 
dynamical  changes  during  seizure  occurrence.  The  overall  goal  of  this  research  is  to  advance  a 
finite  element  based  photoacoustic  tomography  method  for  epilepsy  imaging,  using  both 
laboratory  and  in  vivo  experiments.  Specifically,  in  this  project  we  propose:  (1)  To  design, 
construct  and  test  a  transducer  array  system  for  both  2D  and  3D  PAT  imaging;  (2)  To  advance 
reconstruction  algorithms  and  associated  image  enhancement  schemes  for  quantitative  PAT;  (3) 
To  evaluate  and  optimize  the  integrated  functioning  of  the  hardware  and  software  components 
of  the  transducer  array-based  system,  using  simulation  and  phantom  experiments;  (4)  To  test 
and  validate  the  PAT  system  using  a  well  established  animal  model  of  temporal  lobe  epilepsy. 

Body 

This  report  describes  work  accomplished  in  Year  2  (Months  12-24)  of  the  project.  As  outlined  in 
the  approved  Statement  of  Work  (SOW),  the  tasks  during  this  period  of  time  include:  Task  1 
(Months  0-24):  Assemble  the  entire  photoacoustic  tomography  (PAT)  system;  Task  2  (Months 
18-24):  Design,  build  and  test  the  animal  interfaces;  thus  the  PAT  system  is  ready  for  in  vivo 
studies;  Task  3  (Months  12-24):  Implement  reconstruction  codes  and  associated  enhancing 
schemes  including  dynamic  dual  meshing,  initial  parameter  optimization  and  reconstruction  of 
absolute  optical  absorption  coefficient  for  quantitative  high  resolution  functional  PAT;  and  Task 
4  (Months  12-24):  Calibrate  the  imaging  system;  conduct  simulation  and  extensive  phantom 
experiments  using  the  system  for  evaluation  of  the  reconstruction  codes  and  associated 
enhancements. 

The  sections  below  consist  of  (1)  system  development  and  evaluation  and  (2)  software 
development  that  reflect  the  tasks  associated  with  the  SOW  during  Months  12-24. 


1.  System  Development  and  Evaluation  (Tasks  1,  2  and  4) 

We  are  pleased  to  report  that  we  have  successfully  constructed,  tested  and  calibrated  the 
proposed  transducer  array  based  real-time  PAT  system  in  Year  2.  Here  we  provide  detailed 
description  of  this  novel  PAT  system,  and  present  the  calibration  and  experimental  studies  using 
this  system.  Figs,  la  and  1b  show  the  photograph  of  the  completed  PAT  system. 
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Figure  1.  Photograph  of  the  real-time  PAT  system,  (a):  Top  view,  (b):  Front  view. 


As  shown  in  Fig.  la,  a  tunable  Ti:Sapphire  laser  (part  A)  is  used  to  provide  690-950nm  near- 
infrared  (NIR)  laser  pulses  with  a  full-width  half-maximum  (FWHM)  value  of  8-25ns  pulse  width 
with  a  repetition  rate  of  10Hz  for  generating  photoacoustic  signals.  The  5  MHz  full-ring 
transducer  array  (part  B)  is  used  for  recording  the  photoacoustic  signals.  A  home  made  pre¬ 
amplifier  (part  C)  receives  the  signals  from  the  transducers  and  transmitted  the  amplified  signals 
to  A/D  board  (part  D)  which  are  then  stored  in  the  computer  (part  E).  The  three  key 
components  of  the  system  including  the  tunable  laser,  ultrasound  transducer  array  and 
acquisition  system  are  discussed  in  detail  in  Sections  1.1-1 .3  below. 

Figure  2  depicts  the  system  block  diagram  which  allows  better  understanding  of  the  data 
acquisition  flow  in  this  PAT  system.  The  Ti:Sapphire  laser  optically  pumped  by  a  Q-switched 
Nd:YAG  laser  delivers  8  to  25-ns  pulses  at  10  Hz  with  a  wavelength  tunable  from  690  to  950  nm 
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to  generate  photoacoustic  signals.  The  192-element  transducer  array  is  used  to  capture  the 
photoacoustic  signals.  In  the  full  192-element  configuration,  16  preamplifier  boards,  each 
supporting  12  channels,  amplify  the  signals  from  the  transducer  elements  with  programmable 
gains  of  26  dB.  The  multiplexed  signal  is  subsequently  amplified  with  an  AD604  amplifier  with 
programmable  gain  of  0  to  80  dB,  and  then  digitized  with  the  8x8  channels  A/D  board.  The 
timing  for  data  acquisition  is  specifically  shown  in  Figure  3.  We  see  that  the  timer/controller 
contains  a  10  Hz  laser  Q-switch  signal  repeated  each  100ms  (the  Q-switch  pulse  width  is  lOOus, 
and  laser  pulse  width  is  10ns).  Three  scan  enabling  signals,  ‘EN_CH1-64’,  ‘EN_CH65-128’  and 
‘EN_CH129-192’  are  synchronized  by  the  Q-switch  signal  with  a  3.33Hz  frequency  (period 
300ms).  EN_CH1-64  enables  the  #1~#64  channels  pre-amplifier  and  data  acquisition.  Each 
channel  has  a  sample  rate  of  50MHz,  controlled  by  signal  ‘CLK50M’  (frequency  variable).  The 
post  delay  and  the  sample  waveform  length  are  variable  and  can  be  assigned  by  Labview 
program.  EN_CH65-128  enables  the  #65~#128  channels  pre-amplifier  and  data  acquisitions. 
EN_CH129-192  enables  the  #129~#192  channels  pre-amplifier  and  data  acquisitions.  This 
provides  the  opportunity  to  read  all  captured  pointers  into  the  computer  through  the  local  bus 
after  the  acquisition  and  later  easy  readout  of  the  photoacoustic  data  from  all  channels. 


Ti:  Sapphire 
laser 


Data  Acquisitioi  v_ 


Control  &  PC 


1 6  pre-Amp  boards  in  the  box 


-  1-192 
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Figure  2.  Block  diagram  of  the  real-time  PAT  system. 
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Figure  3  Diagram  of  timing  in  the  PAT  system. 


1.1  Laser  Source 

A  Ti:Sapphire  (Symphotics  Til,  LT-2211A)  laser  optically  pumped  with  a  Q-switched  Nd:YAG 
laser  (Symphotics-TII,  LS-2137/2)  delivers  8-30  ns  pulses  at  10  Hz  with  wavelength  tunable 
from  690  to  1000  nm.  The  beam  is  expanded  with  a  Galilean  telescope  assembly  and 
subsequently  diverged  with  either  a  plano-concave  lens  or  homogenized  by  a  circular  profile 
engineered  diffuser  (ED1-S20,  ThorLabs,  Newton,  NJ).  The  laser  beam  is  positioned  at  the 
center  of  the  transducer  and  strikes  the  stage  orthogonal  to  the  imaging  plane  of  the  transducer 
for  maximum  uniformity. 

1.2  Full-ring  Ultrasound  Transducer  Array 

Based  on  the  extensive  simulation  and  phantom  experiments  in  Year  1,  we  eventually  built  the 
transducer  array,  custom  manufactured  by  Blatek,  Inc.  (PA,  USA),  consists  of  192  elements 
arranged  in  a  full  circular  path  (360  degrees)  with  a  radius  of  150  mm  (Figure  4).  The  transducer 
array  has  a  central  frequency  of  5  MHz  and  a  bandwidth  greater  than  75%.  Each  element  is 
spaced  with  a  lateral  pitch  of  8  mm  and  has  an  active  area  of  3  mm  diameter. 


(a)  Full-ring  transducer  array  (c)  Ultrasound  signals 

Figure  4  Full-ring  ultrasound  transducer  array. 
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1.3  Control  Electronics  and  Data  Acquisition 

The  192  channels  photoacoustic  data  acquisition  system  is  composed  of  4  pre¬ 
amplifier/multiplexer  boxes  and  8  PCIAD850  data  acquisition  boards  (Figure  5).  Each  pre¬ 
amplifier/multiplexer  box  includes  4  boards  with  4x12  AD8099  chips  and  4x4  MAX4051  chips. 
The  whole  system  has  192  AD8099  (26dB)  and  64  MAX4051  chips.  Each  PCIAD850  has  8 
channels  with  a  sample  clock  up  to  50MHz  and  1 0bit  A/D  converter.  The  system  can  sample  64 
channels  at  one  time.  The  #1  PCIAD850  synchronizes  with  the  laser  Q-switch  sync-signal.  The 
computer  sends  out  the  scan  enabling  signals  by  the  USB-1 0241s. 


Figure  5  Flow  diagram  of  the  192  channels  photoacoustic  data  acquisition  system. 

1.3.1  The  first  order  26dB  pre-amplifier 
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In  the  first  order  26dB  pre-amplifier,  the  AD8099  was  choose  to  amplify  the  signal  (Fig.  6a).  In 
the  circuit  path,  R1  should  be  50ohm,  Rf  499ohm,  Rg  26ohm  (Fig.  6b). 


i  ^12  transducers 
/pre-ampNfEer 


4  channel 
Output  after 
pre-amplifier 


(A)  One  of  the  16  p re -amplifiers 


(B)  12  channel 
pre-amplifiers  (2Gdb) 


Figure  6.  (a)  Photograph  of  the  26dB  pre-amplifier,  (b)  26dB  AD8099  design. 


1.3.2  The  second  order  26dB  pre-amplifier  and  52dB  pre-amplifier 

As  a  backup  plan  in  case  of  an  ultra-low  input  signal,  a  52dB  gain  preamplifier  has  been 
designed  (Figure  7). 


1.3.3  Multiplexer  MAX4051 

The  MAX4051  is  used  as  a  3  tol  switch  having  a  ±5  V  DC  power  supply  (Figure  8A-E).  The 
addition  of  diodes  (Figure  9E)  reduces  the  analog  signal  range  to  one  diode  drop  below  V+  and 
one  diode  drop  above  V-,  without  affecting  the  device’s  low  switch  resistance  and  low  leakage 
characteristics.  Here  the  difference  between  V+  and  V-  should  not  exceed  17 V. 
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Figure  8  Design  of  the  MAX4051  as  a  3  tol  switch. 
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1.3.4  The  pre-amplifier/multiplexer  board 

The  whole  system  has  4  pre-amplifier  and  multiplexer  boxes;  each  box  has  4  pre-amplifier  and 
multiplexer  boards.  On  each  board,  there  are  12  analog  input  channels  and  4  multiplexers 
providing  4  amplified  outputs  (Figure  9).  Each  analog  input  channel  can  have  a  gain  of  26dB  or 
52dB. 


(A)  Mapping  between  Amplifier  channels 
&  multiplexer  channel(4/1 2)  i 

i 

i 


(C)  Connector  Pin  Configuration 


(D)  Circuit  diagram 


UF 


I  UNI  V  r.RS  IT  Y  of 

FLORIDA 


Figure  9.  Schematic  of  the  pre-amplifier  and  multiplexer  PCB  board. 


1.3.5  PCIAD850 

PCIAD850  is  an  analog-to-digital  data  converter  (two  boards)  with  simultaneous  acquisition 
ability.  Each  PCIAD850  supports  8  channels,  see  Figure  10.  At  each  trigger  event  (software 
trigger  or  external  trigger),  all  the  channels  will  simultaneously  convert  the  analog  signals  to 
digital  data  with  user  selected  post  trigger  delay  and  waveform  length.  Other  programmable 
parameters  include  sampling  rate,  trigger  source,  trigger  rate,  gain,  DC  offset,  low  and  high 
pass  filters.  Both  boards  feature  user-selectable  10-bit  resolution  and  up  to  50  megabyte 
samples  per  second.  With  one  board  set  as  the  master  and  remaining  set  to  slave  mode,  all 
channels  on  the  slave  boards  will  start  taking  data  upon  receiving  a  trigger  signal  from  the 
master  board.  The  jitter  between  channels  is  less  than  2  nanoseconds.  PCIAD850  has  low  pass 
filters  with  different  cutoff  frequency:  All  frequency,  16  MFIz,  7.3  MHz,  and  4.8  MHz.  It  also  has 
high  pass  filters  with  different  cutoff  frequency:  4.8  MHz,  1.6  MHz,  0.6  MHz,  and  0.016  MHz. 

The  channel  based  parameters: 
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Channel  1 


Block  Diagram  of  PCIAD850/PCIAD1 650  Board 

Channel  2  Channel  16 
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Figure  10.  Block  diagram  of  PCIAD850/1650 


1.3.6  Labview  Data  Acquisition  Software 

Labview  programming  is  used  to  control  the  entire  data  acquisition.  The  software  is  written  to 
allow  maximum  flexibility  for  various  system  configurations.  Pointers  created  during  the  data 
acquisition  process  are  used  to  quickly  access  the  photoacoustic  data.  Data  is  then  saved  for 
image  reconstructions.  Figure  1 1  shows  the  front  control  panel  of  the  Labview  software. 
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Figure  1 1 .  Front  control  panel  of  the  Labview  software. 

1.4  System  Evaluation  and  Experimental  Studies 

1.4.1  System  Calibration 

The  system  calibration  is  primarily  concerned  with  the  calibration  of  the  positioning  and 
response  of  the  multiple  transducers  in  the  system.  We  used  a  photoacoustic  (PA)  point  target 
centrally  located  to  achieve  the  goal.  In  the  calibration,  the  PA  signal  was  measured  by  each 
transducer  at  all  the  locations  along  the  circular  path.  Figure  12  illustrates  the  measured  PA 
signals  from  6  representative  transducers  over  time-of-flight  before  and  after  the  calibration. 
Before  calibration,  we  see  that  the  signal  from  the  point  target  arrives  each  transducer  at 
different  time  indicating  positioning  error  (blue  curves  in  Fig.  12).  By  adjusting  the  position  of 
each  transducer,  we  obtained  the  signals  from  the  point  target  that  arrived  at  the  transducer  at 
the  same  time  point  (red  curves  in  Fig.  12).  Figure  13  shows  the  reconstructed  image  of  the 
three  circular  targets  before  and  after  the  calibration.  We  see  that  the  calibration  provides  much 
improved  image  reconstruction  inters  of  the  target  shape  and  size. 

1.4.2  Spatial  Resolution 

To  evaluate  the  system  spatial  resolution,  we  constructed  a  phantom  consisting  of  two  point 
targets  (brass  wires;  100  pm  in  diameter  each)  suspended  in  tissue  phantom  (solidified  Intralipid 
scattering  medium).  Figure  14a  shows  the  reconstructed  image  of  the  two  point  targets,  while 
Figure  14b  is  the  normalized  profile  of  the  reconstructed  image  along  a  line  y=10  mm.  From  the 
profile,  the  half-width-at-half-maximum  (AB)  and  the  half-width-at-quarter-maximum  (CD)  were 
calculated  to  be  0.13mm.  Using  the  resolution  criterion,  the  spatial  resolution  of  the  system  was 
estimated  to  be  AB+CD-target  diameter=0. 13+0. 13-0. 1=0. 16mm. 
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Before  calibration 
After  calibration 


time  f  tus 

Figure  12.  Measured  PA  signals  from  6  representative  transducers  before  and  after  the  calibration. 


153.460000  mm  151.490000  mm 


Before  calibration 


After  calibration 


Figure  13  Reconstructed  photoacoustic  images  of  three  circular  targets  before  and  after  calibration. 


151.530000  mm 


151.500000  mm.  AB=0. 130000  mm.  CD=0.130000  mm,  resolution=0.160000  mm 


Figure  14  (a)  Reconstructed  image  of  two  point  targets.  (B)  The  normalized  profile  of  the  reconstructed  image  shown 
in  (A)  along  y=10  mm. 
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1.4.3  Phantom  Evaluation 

Phantom  experiments  that  simulate  detection  of  seizure  focus  have  been  conducted.  Fig.  15a 
shows  the  recovered  photoacoustic  image  of  O.lmm-diameter  target  mimicking  a  single  seizure 
focus,  while  Fog.  15b  gives  the  reconstructed  image  of  three  targets  having  a  diameter  of  4  mm, 
1.2  mm  and  0.6  mm,  simulating  multiple  seizure  foci.  All  the  images  obtained  show  that  these 
targets  with  different  sizes  can  be  clearly  detected  using  our  imaging  system. 


Figure  15  Reconstructed  image  of  a  single  target  (a)  and  three  targets  (b). 

Figure  16  presents  cross-sectional  photoacoustic  images  of  a  simulated  blood  vasculature 
located  at  various  depths  ranging  from  0  to  4.6  cm  below  the  phantom  surface.  The  cross- 
sections  illustrate  the  capability  of  the  system  to  provide  high  penetration  imaging  depth  with 
high  resolution  and  contrast  at  least  at  a  depth  up  to  3.5  cm  (Figures  16a  and  16b  show  the 
photograph  of  the  samples). 


Figure  16  Cross-sectional  images  of  simulated  blood  vasculature  at  various  depths. 


We  also  conducted  phantom  experiments  to  test  the  high  temporal  resolution  ability  of  the 
system  for  real-time  tomographic  imaging  of  dynamic  fluid  flow  (3-1 0Hz).  Figure  17  depicts 
frames  from  a  10-second  sequence  of  dynamic  ink  flow  through  a  1.6  mm-diameter 
polyethylene  tube  using  a  manual  syringe  push.  To  avoid  partial  view  effects  and  maximize 
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temporal  resolution,  a  sliding  window  acquisition  was  adopted  in  which  element  data  from  the 
nearest  (or  equal)  time  instant  were  used  for  image  formation.  The  images  obtained  show  that 
we  can  clearly  track  the  flow  through  the  tubing  with  high  spatial  and  temporal  resolution.  Such 
tomographic  imaging  dynamic  fluid  flow  was  previously  impossible.  This  will  provide  a  powerful 
tool  for  in  vivo  real-time  monitoring  of  brain  hemodynamics  during  epileptic  seizure  onset. 


laser 


Tube 


A.  Experiment  diagram 


C.  Real  time  PA  movie 


8.  PA  slices  during  circulation 

I  TT^  I  UNIVERSITY  of 

UF I  FLORIDA 


Figure  17  Photoacoustic  images  of  dynamic  ink  flow  through  a  1 .6  mm  diameter  tube. 


1.4.4  Animal  Interface  and  Initial  In  Vivo  Testing  of  the  System 

We  have  constructed  the  animal  interface,  an  animal  chamber  with  adjustable  vertical  position 
with  respective  to  the  transducer  array  (see  Fig.  18).  During  an  experiment,  the  rat  head  is 
placed  into  the  imaging  region  through  an  opening  covered  with  a  piece  of  polyethylene 
membrane  at  the  top  of  the  animal  chamber. 

We  have  performed  initial  animal/rat  experiments  to  test  the  effectiveness  of  the  animal 
interface  for  in  vivo  imaging.  Figure  19a  shows  the  cross-sectional  photoacoustic  image  of  a 
normal  rat  brain  using  the  interface  and  the  PAT  system,  while  Figure  19b  gives  the  open  skull 
photograph  of  the  same  rat  after  the  PAT  imaging  was  performed.  Compared  to  the  open  skull 
photograph,  we  see  that  major  structures  and  various  vasculature  of  the  rat  brain  are  clearly 
imaged  with  high  quality.  Figure  20  presents  the  recovered  photoacoustic  images  of  another  rat 
brain  at  various  depths  ranging  from  1  and  9  mm  below  the  skin  along  with  the  open-skull 
photograph  viewed  from  dorsal  and  basal  surfaces.  At  the  dorsal  plane,  the  superior  sagittal  and 
cerebral  veins  including  branches  are  visible  in  addition  to  the  eyes.  At  the  deeper  cross- 
sections,  the  basalar  veins  and  the  middle  cerebral  artery  can  be  observed  along  with  ringed 
features  near  the  brain  stalk  not  visible  in  the  photograph.  This  demonstrates  that  our 
photoacoustic  imaging  system  has  the  potential  to  provide  3D  imaging  of  the  rat  brain  for 
accurate  epileptic  seizure  localization. 
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Figure  18  The  photograph  of  the  animal  interface  and  transducer  array. 


Min 


Optical  absorption 


Max 


(B)  Skull  opened 


(A)  Noninvasive  PA  brain  image  in  vivo 

Figure  19  (a)  Recovered  photoacoustic  image  for  a  normal  rat  brain.  The  color  scale  (below) 
represents  the  optical  absorption  of  tissue  (arbitrary  units),  (b)  Photograph  of  the  rat  cerebral  with 
the  skull  opened. 
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(C)  Cross-sectional  images  of  rat  brain  vasculature  at  various  depths 

Figure  20  Recovered  photoacoustic  images  of  a  normal  rat  brain  at  various  depths  below  the  skin. 


2.  Software  Development  (Tasks  3  and  4) 

In  Task  3,  we  proposed  to  implement  6  image  enhancement  schemes.  In  Year  1,  3  schemes 
including  total  variation  (TV),  advanced  regularization  techniques,  and  adaptive  meshing  were 
developed  successfully.  In  Year  2,  the  remaining  3  schemes,  including  dynamic  dual  meshing, 
initial  parameter  optimization,  and  reconstruction  of  absolute  optical  absorption  coefficient  for 
quantitative  high  resolution  functional  PAT  have  been  implemented  as  proposed. 

2.1  Dynamic  dual  meshing 

The  concept  of  dual  meshing  is  to  decouple  the  spatial  parameterization  of  the  image  property 
to  be  reconstructed  from  the  model-based  solution  for  the  measured  quantity.  This  allows 
specification  or  discretization  of  the  reconstruction  parameter  to  be  independent  of  the  model 
solution  which  is  a  powerful  approach  in  that  it  lets  the  sampling  density  of  these  two  quantities 
be  distinct.  Hence,  the  use  of  dual  mesh  method  can  considerably  reduce  computational  costs 
during  image  formation  while  improving  the  accuracy  of  image  reconstruction. 

In  Year  2,  we  have  developed  the  proposed  dynamic  dual  meshing  scheme  that  is  able  to 
automatically  solve  the  problem  of  mesh  overlapping  encountered  for  epilepsy  imaging.  In 
particular,  the  new  scheme  can  achieve  adaptive  dual  meshing  that  uses  information  about  the 
resultant  image  as  it  evolves  during  the  iterative  reconstruction  process.  Consequently  it  has  the 
capability  to  adapt  the  dual  mesh  node  deployment  into  a  more  optimal  spacing  for  the  specific 
image  under  construction. 

Two  phantom  experiments  were  conducted  to  test  the  dynamic  dual  meshing  scheme.  We 
employed  single-target-containing  phantom  tests,  aiming  to  validate  the  accuracy  improvement 
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in  detecting  target  when  the  dynamic  dual  mesh  is  used.  The  phantom  materials  used  consisted 
of  Intralipid  as  scatter  and  India  ink  as  absorber  with  Agar  powder  (1-2%)  for  solidifying  the 
Intralipid  and  India  ink  solution.  The  absorption  coefficient  of  the  background  phantom  was  0.01 
mm'1,  while  the  absorption  coefficient  of  the  target  was  0.015  mm'1  and  0.02  mm'1  for  test  1  and 
test  2,  respectively. 

Figure  21  shows  the  evolution  of  finite  element  mesh  during  inverse  iteration  procedure  when 
the  new  scheme  is  used  for  PAT.  Figure  22  shows  the  reconstructed  photoacoustic  images  for 
the  two  cases  with  uniform  mesh  and  dynamic  dual  mesh.  We  see  that  the  target  is  clearly 
better  recovered  with  the  dynamic  dual  mesh  (a  &  d)  in  terms  of  its  shape  and  size  over  that 
with  the  uniform  mesh  (c  &  f)  for  both  cases.  We  also  note  that  the  background  region  for  both 
cases  is  more  smoothly  recovered  using  the  new  scheme  compared  to  that  using  the  uniform 
mesh. 


(a)  (b)  (c) 


Fig.  21  (a)  The  initial  mesh  used  for  image  reconstruction;  (b)  Dynamic  adaptive  mesh  generated  after  the  1st  iteration 
with  one-time  refinement;  (c)  Dynamic  mesh  generated  after  the  2nd  iteration. 


(a)  (b)  (c) 


(d)  (e)  (f) 

Fig.  22  Reconstructed  photoacoustic  images  from  experimental  data,  (a):  the  image  for  case  1  with 
dynamic  dual  mesh,  (b):  the  image  for  case  1  with  dynamic  dual  mesh  after  1  iteration,  (c):  the  image  for 
case  1  with  uniform  mesh,  (d):  the  image  for  case  2  with  dynamic  dual  mesh,  (e):  the  image  for  case  2 
with  dynamic  dual  mesh  after  1  iteration,  (c):  the  image  for  case  2  with  uniform  mesh. 


2.2  initial  parameter  optimization 


16 


Huabei  Jiang,  PhD/Annual  Report 


It  has  been  shown  in  optical  image  reconstruction  that  we  have  realized  a  novel  scheme  to 
optimize  the  boundary  conditions  (BCs)  coefficients  as  well  as  other  initial  parameters  including 
the  initial  estimates  of  the  optical  properties.  Based  on  this  experience,  in  Year  2,  we  have 
implemented  a  new  scheme  to  optimize  the  initial  parameters  needed  for  image  reconstruction 
in  PAT,  which  utilizes  the  parallel  structure  that  is  able  to  significantly  reduce  the  computational 
costs  for  the  data  pre-processing. 

Fig.  23  shows  a  set  of  typical  reconstructions  for  the  optimization  of  initial  parameters  study 
where  the  acoustic  velocity  is  varied  from  1.0x106to  1.6*106  mm/s  for  a  test  case  having  one 
target  embedded  in  a  circular  background  medium.  We  note  that  the  optimized  initial 
parameters,  particularly  the  acoustic  velocity  is  useful  in  reducing  the  noise  effect  in  the 
background  and  improve  the  quantitative  nature  of  the  recovered  absorbed  energy  density. 


(a)  (b)  (c)  (d) 

Fig.  23  (a):  Exact  distribution  map  of  absorbed  energy  density,  (b):  Recovered  energy  density  map  with  velocity 
equal  to  1 .0x1 ,06  (b),  1 .4x1 ,06  (c),  and  1 .6x1 ,06  mm/s  (d). 

2.3  Reconstruction  of  absolute  optical  absorption  coefficient  for  quantitative  high  resolution 
functional  PAT 

One  major  limitation  of  the  existing  method  for  reconstructing  absorption  coefficient  is  that  it 
requires  the  measurement  of  the  absolute  excitation  light  strength.  In  addition,  it  does  not  use 
regularization  techniques,  thus  the  inverse  solution  procedure  sometimes  is  not  stable.  To 
overcome  these  limitations,  we  have  developed  a  novel  reconstruction  method,  where  images 
of  optical  absorption  coefficient  are  obtained  from  a  diffusion  equation  based  regularized 
Newton  method  with  the  absorbed  energy  density  distribution  from  conventional  photoacoustic 
tomography  serving  as  the  measured  field  data.  To  recover  the  absorption  coefficient  from  the 
absorbed  energy  density,  o> ,  the  photon  diffusion  equation  as  well  as  the  Robin  boundary 
conditions(BCs)  can  be  written  in  consideration  of 

V-D(r)V(E(r)$>(r))-$>(r)  =  -S(r)  (1) 

-DV(E(r)<S>)-n  =  E(r)cc<S>  (2) 

where  /Ja  is  the  optical  absorption  coefficient,  'Fis  the  photon  density,  E(r)  =  \/  jujr) ,  D(r)  is  the 
diffusion  coefficient,  D  =  \/(3(jua  +ju's))  and  ju's  is  the  reduced  scattering  coefficients,  a  is  a 
boundary  condition  coefficient  related  to  the  internal  reflection  at  the  boundary,  and  S(r)  is  the 
incident  point  or  distributed  source  term.  For  the  inverse  computation,  the  initial  estimate  of 
inverse  of  absorption  coefficient  can  be  updated  based  on  iterative  Newton  method  as  follows 

A(E)  =  ( Jr J  +  TI  +  LrL)_1  [ Jr  (<I>°  -<PC)]  (3) 

in  which  L  is  the  filter  matrix,  A  is  the  regularization  parameter  and  J  is  the  Jacobian  matrix.  We 
experimentally  demonstrate  this  new  method  using  tissue-mimicking  phantom  measurements. 
For  the  two  experiments,  the  phantom  materials  used  consisted  of  Intralipid  as  scatterer  and 
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India  ink  as  absorber  with  Agar  powder  (1-2%)  for  solidifying  the  Intralipid  and  India  ink  solution. 
The  background  phantom  had  =0.01  mm'1  and  //j  =  1.0mm'1  while  the  two  targets  had 
//a=0.03  mm'1  and //’=2.0  mm"1  for  test  1,  and  /ua=0.07  mm'1  and  ju's  =3. 0mm"1  for  test  2.  Figs. 
24(a)  and  24(b)  present  the  reconstructed  absorption  coefficient  images  of  two  objects  having  a 
size  of  2.0  and  3.0mm  (test  1),  and  5.5mm  (test  2)  in  diameter,  respectively,  while  the  recovered 
absorbed  energy  density  maps  for  experiments  1  and  2  are  plotted  in  Figs.  24(c)  and  24(d)  for 
comparison.  The  reconstruction  results  show  that  the  optical  absorption  coefficient  images 
obtained  are  quantitative  in  terms  of  the  shape,  size,  location  and  optical  property  values  of  the 
heterogeneities  examined. 


(a) 


I 


|0.O7 

.06 

■0.05 

0.04 

003 

0,02 


(C) 


0.02 

0.015 

a  oi 

0.005 


(d) 


Fig.  24  Reconstructed  absorption  coefficient  images  (a,  b)  and  absorbed  optical  energy 
density  images  (c,  d)  for  tests  1  and  2.  (a)  and  (c)  are  for  test  1 ,  and  (b)  and  (d)  for  test  2. 


Key  Research  Accomplishments 

1 .  We  have  completed  the  construction  of  the  proposed  an  array  based  real-time  PAT  system, 
and  calibrated  and  tested  the  system  using  extensive  phantom  experiments. 

2.  We  have  built  the  animal  interface  and  successfully  tested  it  for  in  vivo  imaging  of  rat  brain. 

3.  We  have  developed  three  novel  schemes  that  can  enhance  our  current  reconstruction 
software. 

4.  We  have  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 


Reportable  Outcomes  (see  the  Appendix  to  this  Summary  Report) 


We  expect  that  the  PAT  system  completed  in  Year  2  will  allow  us  to  generate  considerable  data 
for  writing  numerous  journal  publications  and  conference  proceedings  in  Year  3. 
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Conclusions 

We  have  made  a  significant  progress  that  has  fulfilled  the  statement  of  work  proposed  for  Year 
2  of  this  project.  The  real-time  PAT  system  completed  in  Years  1  and  2  provides  us  a  platform 
for  performing  the  extensive  phantom  and  animal  experiments  proposed  for  the  coming  years. 
We  are  confident  that  we  will  be  able  to  fulfill  or  exceed  the  work  statement  for  Years  3  and  4. 

Appendix 

L.  Yao  and  H.  Jiang,  Enhancing  finite  element  based  photoacoustic  tomography  using  total- 
variation  minimization,  Physics  in  Medicine  and  Biology  (submitted). 
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Abstract 


A  total-variation-minimization-based  finite  element  reconstruction  algorithm  for 
photoacoustic  tomography  is  described  in  th  is  paper  that  enhances  the  quality  of 
reconstructed  photoacoustic  im ages  with  tim e-domain  data.  Simulations  are  conducted 
where  different  contrast  levels  betw  een  the  target  and  the  bac  kground,  different  noise 
levels,  different  sizes  of  and  shapes  of  the  ta  rget  are  studied  in  a  30mm-diameter  circular 
heterogeneous  background.  These  sim  ulated  results  show  that  the  quality  of  the 
reconstructed  images  can  be  im  proved  significantly  due  to  the  dec  reased  sensitivity  to 
noise  effect  when  the  tota  1-variation  minimization  is  included.  The  enhancem  ent  is 
further  confirmed  using  experim  ental  data  obtained  from  several  phantom  experiments 
and  an  in  vivo  animal  experiment. 
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1.  Introduction 


Photoacoustic  tomography  (PAT)  is  an  em  erging  non-invasive  im  aging  technique  that 
combines  the  merits  of  high  optical  contra  st  and  high  ultrasound  resolution  in  a  single 
modality  [1-3].  In  PAT,  the  Helm  holtz-like  photoacoustic  (PA)  wave  equation  has  been 
commonly  used  as  an  accurate  model  for  describing  laser-induced  acoustic  wave 
propagation  in  tissue.  While  analytical  re  construction  methods  have  been  used  for 
photoacoustic  image  reconstructions,  the  fin  ite  element  method  (FEM)  based  approach 
appears  to  be  particularly  powerful  in  this  regard  [4-5].  The  advantages  of  the  FEM  based 
PAT  method  include  (1)  quantitative  imaging  capability  by  recovering  optical  absorption 
coefficient,  (2)  elimination  of  the  assumption  of  homogeneous  acoustic  medium  needed 
in  analytical  m  ethods,  (3)  accommodation  of  object  boundary  i  rregularity,  and  (4) 
appropriate  boundary  conditions  implementations. 

Our  finite  elem  ent  PAT  approaches  are  b  ased  on  a  regularized  iterative  Newton 
method,  which  have  been  tested  and  evaluate  d  using  extensive  simulations  and  phantom 
experiments  in  both  f  requency-  and  tim  e-domain  [4-7].  While  these  re  suits  are 
encouraging  and  promising,  we  realize  that  m easurement  noises  are  still  the  major  factor 
affecting  the  quantitative  accuracy  of  the  recon  structed  images.  For  example,  erro  rs  for 
quantitatively  recovering  optical  absorption  coefficient  can  be  as  large  as  10%  for 
simulated  data  with  5%  noise  added  and  20%  for  experimental  data  [6-7]. 

In  an  effort  to  reduce  the  noise  effect  and  to  further  enhance  the  quantitative  accuracy 
of  photoacoustic  im  age  reconstruction,  in  th  is  paper  we  consider  a  total-variation- 
minimization  scheme  within  our  FEM-base  d  reconstruction  fram  ework.  Our  existing 
reconstruction  algorithms  are  based  on  the  leas  t-squares  criteria  (i.e.,  the  regula  rized 
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Newton  method)  [11-12]  that  stand  on  the  sta  tistical  argument  that  the  least-squares 
estimation  is  the  bes  t  estimator  over  an  entire  ensem  ble  of  all  possib  le  pictures.  Total 
variation,  on  the  other  hand,  m  easures  the  oscillations  of  a  given  function  and  does  not 
unduly  punish  discontinuities  [8-9].  Hence,  one  can  hypothesize  that  a  hybrid  of  these 
two  minimization  schemes  should  be  ab  le  to  provide  higher-quality  image 
reconstructions.  In  fact,  the  st  rategy  of  finding  m  inimal  total-variation  has  proved  to  be 
successful  in  applications  including  elec  trical-impedance  tomography  [9],  microwave 
imaging  [10],  im  age  processing  [8,  13-14],  optimal  design  [15]  and  diffuse  optica  1 
tomography  [16]. 

In  the  f  ollowing  sections,  we  will  describ  e  in  detail  the  im  plementation  of  total- 
variation-minimization  within  our  existi  ng  PAT  reconstruction  fram  ework  in  tim  e- 
domain.  We  will  dem  onstrate  the  enhanced  reconstruction  algorithm  using  sim  ulated, 
phantom  and  in  vivo  data. 

2.  Total- Variation-Minimization  Scheme 

To  describe  our  total- variation-m  inimization  (TVM)  method,  we  first  briefly  introduce 
the  FEM-based  regularized  Newton  m  ethod  in  time-domain  [7],  The  tim  e-domain 
photoacoustic  wave  equation  in  tissue  can  be  described  as  follows 

(ij 

V  '  v02  8t2  cpdt 

where  p  is  the  pressure  wave;  v0  is  the  speed  of  acoustic  wave  in  the  medium;  /?  is  the 
thermal  expansion  coefficient;  Cp  is  the  specific  heat;  <J>  is  the  absorbed  energy  density; 
j(t)  =  S(t-t0)  is  assumed  in  our  study. 

The  finite-element  discretization  of  Eq.  (1)  can  then  be  written  as 
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2>, 

7=1 


r  A  r  1 

J V^,  -Vi//jdS  +Y,Pj  \s—¥iYjdS 

S  j= 1  L  vo 


where  y/ .  is  the  basis  function. 


The  first-order  absorbing  boundary  conditions  used  are  [17-20]: 

1  dp  p 


Vp-n  = 


v0  dt  2  r 


(3) 


where  h  is  the  unit  normal  vector. 

In  both  the  forward  and  inverse  calculations,  the  unknown  coefficient,  O  needs  to 
be  expanded  in  a  sim  ilar  fashion  to  p  as  a  sum  of  unknown  param  eters  multiplied  by  a 
set  of  locally  spatially  varying  Lagrange-p  olynomial  basis  functions.  Thus  the  matrix 
form  of  Eq.  (2)  becomes, 

MMdcM+Mp}=M>  (4) 

where  the  elements  of  matrix  [/f],  [c]  and  [m]  are 


Kij=\^ Vi  ■  V  VjdS  +  ^§iy/,y/Jdl ; 

S’  ^ V 

cij  =^iv,Vjdl; 

Mii=\\sV,VjdS ; 


and  the  column  vectors  \p\,  [p],  [p]  and  {5}  are 


c 


dS-  —  ; 
dt 


{p)={Pi,P2,---Pn}T ;  {p}={PmP2>---Pn}T  {p}={Pi,P2>---Pn}T ’ 
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Here  the  Newmark’s  time-stepping  scheme  has  been  used  for  the  discretization  of  time 
dimension  [21-22],  which  is  a  commonly  used  im  plicit  method  for  the  second-order 
propagation  equations  such  as  Eq.  (4). 

To  form  an  image  from  a  presumably  uniform  initial  guess  of  the  absorbed  en  ergy 
density  distribution  we  need  a  method  of  updating  O  from  its  starting  value.  This  update 
is  accomplished  through  the  least-squares  minimization  of  the  following  functional: 

M  /  X, 

F(p,<l>)=UP°-Pi)  ’  (5) 

j= 1 

where  p’  and  pCj  are  observed  and  computed  acoustic  field  data  for  i  =  1,2,- •  • M 

boundary  locations.  Using  the  regularized  Newton  method,  we  obt  ained  the  following 
equation  for  updating  O  : 

(3r3  +  ^  =  3r(/-P').  («) 

where  p°  =(p°,  P°2,-"P°UJ  ,  p‘  =  (p,',  p\,  is  the  update  vector  for  the 

absorbed  optical  energy  density,  3  is  the  Jacobian  m  atrix  formed  by  dp / <3d>  at  the 
boundary  measurement  sites;  X  is  the  regularization  parameter  determined  by  combined 
Marquardt  and  Tikhonov  regularization  schemes,  and  I  is  the  identity  matrix. 

Two  typical  approaches  ex  ist  for  minimizing  total  va  riation:  a  constrained 
minimization  through  the  solution  of  the  nonlinear  PA  equation  [8,  23]  and  an 
unconstrained  minimization  by  addition  of  the  total  variation  as  a  penalty  term  to  the  least 
squares  functional  [10,  13-14,  24],  Solutio  ns  obtained  with  th  ese  two  m  ethods  are 
essentially  the  sam  e;  however,  from  a  computational  standpoint,  unconstrained 
minimizations  are  much  easier  to  implem  ent  and  require  much  less  modifications  to  the 
existing  algorithm  [13].  Thus  in  this  study  the  unconstrained  TVM  is  used. 
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We  now  incorporate  the  total  variation  of  0  as  penalty  term  by  de  fining  a  new 
functional  [9-10,  13]: 

F{p,®)  =  F{p,®)+L(®).  (7) 

Here 


+  S2  dxdy ,  (8) 

is  the  penalty  term,  and  5  are  typically  positive  parameters  that  need  to  be  determined 
numerically.  The  m  inimization  of  Eq.  (7)  proceeds  in  standard  fashion  by  the 

differentiation  of  F  with  respect  to  each  nodal  para  meter  that  constitutes  th  e  <J> 
distribution;  simultaneously  all  these  relations  are  set  to  zero.  These  steps  lead  to  the 
following  nonlinear  system  of  equations 
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Similarly  to  Eq.  (6),  the  following  matrix  equation  for  TVM  constrained  inversion  can  be 
obtained 

{^TZ  +  R  +  My±x  =  ZT{p° -pc)-V ,  (10) 


where 
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and  V  =  {Vl,V2,-VN}T. 

3.  Results  and  Discussion 

In  this  section  the  TVM  enhanced  reconstructi  on  algorithm  is  tested  and  evaluated  using 
both  simulated  and  experimental  data.  For  comparative  purposes,  reconstructions  without 
the  TVM  enhancement  are  also  presented. 

3.1  Simulations 

In  these  simulations,  a  dual  meshing  method,  as  detailed  elsewhere  [6],  is  used  for  fast 
yet  accurate  inverse  com  putation.  The  fine  mesh  used  for  the  forward  calcu  lation 
consisted  of  3627  nodes  and  7072  e  lements,  while  the  coarse  m  esh  used  for  the  inverse 
calculation  had  930  nodes  and  1768  elements.  All  the  images  obtained  from  the  method 
without  the  TVM  are  the  results  of  three  it  erations,  while  those  obtained  from  the  TVM- 
enhanced  method  are  the  results  of  fifteen  iterations.  (W  e  found  that  larger  num  her  of 
iteration  changed  the  solutions  only  by  less  than  0.5%).  For  the  sim  ulations  presented, 
&><(,  =  0.5  and  8  =  0.001  for  the  first,  second  and  third  cases,  and  co<b  =  1.0  and  8  =  0.001 
for  the  fourth  case  were  used. 

For  the  first  simulation,  the  test  g  eometry  is  a  30mm-diam  eter  circular  background 
containing  an  off-center  4-  mm-diameter  target  region.  The  target  had  0  =  1.0  mJ / mmi  , 
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while  the  background  had  O  =  2.0»?y/ mm3  .  In  this  case  the  image  recons  traction  was 
performed  with  0%,  10%  and  25%  noise  added  to  the  “m  easured”  data.  Fig.  1  gives  two 
sets  of  absorbed  energy  density  im  ages  recovered  using  the  reconstruction  m  ethod 
without  (Figs.  la,c,e)  and  with  the  TVM  enhancement  (Figs.  lb,d,f)  under  the  conditions 
of  different  noise  levels.  As  can  b  e  seen,  enhancement  of  the  recon  straction  by  the 
incorporation  of  the  TVM  is  obvious  over  th  e  method  without  the  TV  M.  To  provide  a 
more  quantitative  assessment  of  these  im  ages,  Fig.  2  is  inc  luded,  in  whic  h  the 
reconstructed  absorbed  energy  density  profiles  are  displaye  d  along  transects  through  the 
target  for  the  images  shown  in  Fig.  1.  We  find  that  the  TVM  enhanced  method  not  only 
can  improve  the  quality  of  the  recovered  images,  but  also  can  decrease  the  sensitivity  of 
the  method  to  noise  effect. 

The  second  simulation  is  intended  to  investig  ate  how  the  target  si  ze  affects  the  TVM 
enhancement.  In  this  cas  e,  no  noise  was  added  to  the  “m  easured”  data,  and  the  diam  eter 
of  the  off-center  target  was  2  mm.  Figs.  3a  and  3b  give  two  sets  of  absorbed  energy 
density  images  recovered  using  the  m  ethod  without  the  TVM  (Fig.  3a)  and  with  the 
enhancement  (Fig.  3b),  while  Figs.  3c  and  3d  present  the  quantitat  ive  profiles  of  the 
absorbed  energy  density  along  transects  through  the  target  for  the  images  shown  in  Figs, 
la, lb  and  Figs.  3a, 3b.  Again,  considerable  improvement  can  be  observed  from  the 
reconstructed  results  when  the  TVM  is  i  nvoked  compared  with  the  m  ethod  without  the 
TVM.  It  is  also  interesting  to  note  that  the  enhancement  for  this  case  is  more  striking  than 
that  for  the  first  one  where  the  background  contained  a  larger  target. 

The  third  simulation  aims  to  see  how  the  contrast  level  of  the  absorbed  energy  density 
between  the  target  and  the  backg  round  impacts  the  TVM  enhancement.  In  this  case,  no 
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noise  was  added  to  the  “measured”  data,  and  the  off-center  target  had  a  diameter  of  4mm 
and  <f>  =  1 .5mj / mm2'  .  Figs.  4a  and  4b  present  two  sets  of  absorbed  energy  density  images 
recovered  using  the  method  without  the  TYM  (Fig.  4a)  and  with  the  T  VM  enhancement 
(Fig.  4b),  while  Figs.  4c  and  4d  show  the  qua  ntitative  profiles  of  the  absorbed  en  ergy 
density  along  transects  through  the  target  for  the  im  ages  shown  in  Figs,  la, lb  and  F  igs. 
4a, 4b.  We  can  see  that  the  im  ages  formed  by  incorporation  of  the  TVM  are  clearly 
enhanced  qualitatively  in  visual  content  relative  to  that  obt  ained  using  the  m  ethod 
without  the  TVM.  We  also  note  that  lower  the  contrast  level  is,  larger  the  enhancement  is. 

In  the  fourth  simulation,  three  targets  having  different  shapes  (circular,  elliptical  and 
rectangular)  were  embedded  in  the  background,  and  the  absorbed  energy  density  of  these 
targets  was  2.5  mj/mm*  ,  1 .5  mj /mm3  and  2.0 mj j mm3  ,  respectively.  A  noise  level  of 
25%  was  added  to  the  “m  easured”  data  in  th  is  case.  Fig.  5  shows  the  exact  and  the 
recovered  absorbed  energy  density  im  ages  as  well  as  th  e  quantitative  absorbed  en  ergy 
density  property  profdes  al  ong  the  transect  that  across  these  targets.  Again,  the 
improvement  in  image  quality  is  apparent. 

3.2  Phantom  and  In  Vivo  Experiments 

In  this  section  both  phantom  and  in  vivo  expe  rimental  data  were  used  to  confirm  our 
findings  from  the  si  mulations.  The  experimental  setup  used  for  collecting  the  phantom 
data  was  a  pulsed  ND:  YAG  lase  r  based  single  transducer  (1MHz)  scanning  system, 
which  was  described  in  detail  elsewher  e  [5],  Thre  e  phantom  experim  ents  were 
conducted.  In  the  first  two  experiments,  we  embedded  one  or  two  objects  with  a  size 
ranging  from  3  to  0.5  mm  in  a  50  mm-diameter  solid  cylindrical  phantom.  The  phantom 
materials  used  consisted  of  Intralipid  as  s  catterer  and  India  ink  as  abs  orber  with  Agar 
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powder  (1-2%)  for  solidifying  the  Intralip  id  and  India  ink  solution.  The  absorption 
coefficient  of  the  background  phantom  was  0.01  mm'1,  while  the  absorption  coefficient  of 
the  target(s)  was  0.03  mm  In  the  last  experim  ent,  we  used  a  single -target-containing 
phantom,  aiming  to  test  the  capability  of  de  tecting  target  having  low  optical  contrasts 
relative  to  the  background  phantom.  In  this  case,  the  target  had  an  absorption  coefficient 
of  0.015  mm  The  reduced  scattering  coefficients  of  the  background  phantom  and 
targets  were  1.0  and  3.0  mm'1  for  the  first  two  experiments,  and  1.0  and  2.0  mm'1  for  the 
last  experiment. 

A  total  of  120  receivers  were  equally  dist  ributed  along  the  surface  of  the  circular 
background  region.  In  the  reconstructions,  the  fine  mesh  used  for  the  forward  calculation 
consisted  of  5977  nodes  and  1 1712  elem  ents,  while  the  coarse  mesh  used  for  the  inverse 
calculation  had  1525  nodes  and  2928  elements.  The  reconstructed  images  were  the  results 
of  three  and  fifteen  iterations  for  the  method  without  and  with  the  TVM,  respectively.  For 
these  experimental  cases,  co#  =  1.0  and  8  =  0.001  for  the  first  and  third  cases  and 
atj,  =  2.0  and  8  =  0.001  for  the  second  case  were  used. 

Figure  6  shows  the  reconstruc  ted  absorption  coefficient  im  ages  for  all  the  th  ree 
experimental  cases,  while  Figure  7  presents  quantitative  absorption  coefficient  profiles 
along  transects  through  one  ta  rget  for  the  im  ages  shown  in  Fig.  6.  W  e  see  that 
considerably  enhanced  images  are  achieve  d  with  th  e  total-variation  minimization, 
especially  when  the  target  is  small  (case  2),  or  when  the  contrast  level  between  the  target 
and  the  background  is  low  (case  3). 

As  a  final  exam  pie,  our  methods  were  tested  using  in  vivo  data  collected  previously 
from  an  a  nimal  (rat)  m  odel  of  epileps  y  [25],  Focal  seizures  were  induced  by 
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micro  injection  of  bicuculline  methiodide  (BMI)  into  the  parietal  neocortex.  The  in  vivo 
were  collected  using  the  same  pulsed  ND:  YAG  laser  based  scanning  PAT  system.  In  the 
reconstruction,  the  fine  mesh  used  for  the  forward  calculation  consisted  of  17713  nodes 
and  34944  elements,  while  the  coarse  m  esh  used  for  the  inverse  calculation  had  4489 
nodes  and  8736  elements.  The  in  vivo  im  ages  obtained  were  the  results  of  two  and  five 
iterations  for  the  m  ethods  without  and  with  the  TVM,  respectiv  ely.  We  found  that 
&><(,  =0.1  and  8  =  0.001  appeared  to  provide  optimal  results. 

Fig.  8  gives  the  recovered  images  for  the  rat  brain  scanned  10  min  after  the  injection  of 
BMI  without  and  with  the  use  of  the  TVM.  The  in  vivo  resu  Its  shown  here  confirm  that 
both  of  the  reconstruction  m  ethods  can  provide  quality  im  ages  (The  arrow  indicates  the 
detected  seizure  focus),  while  no  table  enhancement  in  im  age  quality  can  be  observed 
with  the  TVM  based  method.  For  example,  the  actual  continuous  middle  cerebral  artery 
(indicated  by  the  dashed  circle)  is  shown  discontinued  in  Fig.  8a  (without  the  T  VM), 
while  this  feature  is  clearly  depicted  in  Fig.  8b  (with  the  TVM). 

4.  Conclusions 

We  have  presented  a  tim  e-domain  FEM  based  photoacoustic  im  age  reconstruction 
method  that  incorporates  the  TVM.  The  results  shown  in  this  work  indicate  that  the  TVM 
based  method  is  able  to  offer  clear  enhancement  in  image  reconstruction  over  the  method 
without  the  TVM,  not  only  in  te  rms  of  the  location,  size  and  shape  of  the  target,  but  also 
the  absorbed  energy  density  property/optical  absorption  coefficient  values  themselves.  In 
addition,  the  results  have  shown  that  the  inclusion  of  the  TVM  in  our  recon  struction 
algorithm  is  highly  effective  in  the  presence  of  noisy  data.  Finally,  it  is  important  to  note 
that  the  TV  M  based  method  m  ay  be  ideal  for  image  reconstruction  involving  a  low 
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contrast  level  between  the  ta  rget  and  the  background,  or  sm  all  targets  embedded  in  the 
background. 
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Figure  captions 


Fig.l  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and 
without  the  total-variation-minimization  (TYM)  enhancement  under  different  noise  levels 
(case  1).  (a),  without  the  TVM,  0%  noise,  (b),  with  the  TVM,  0%  noise,  (c),  without  the 
TVM,  10%  noise,  (d),  with  the  TVM,  10%  noise,  (e),  without  the  TVM,  25%  noise,  (f), 
with  the  TVM,  25%  noise. 

Fig.2  Comparison  of  the  exact  and  reconstruc  ted  absorbed  energy  de  nsity  profiles  along 
transect  y=0mm  for  the  im  ages  appearing  in  Fig.l.  (a),  0%  noise,  (b),  10%  noise,  (c), 
25%  noise. 

Fig.3  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and 
without  the  total-variation-minimization  (TVM)  enhancement  with  different  target  sizes 
(case  2).  (a),  without  the  TVM,  2mm  -diameter  target,  (b),  with  the  TVM,  2mm  -diameter 
target,  (c),  the  absorbed  energy  density  profiles  along  the  transect  y=0mm  for  the  i  mages 
appearing  in  Figs,  la  and  lb  (4mm  -diameter  target),  (d),  the  absorb  ed  energy  density 
property  profiles  along  the  transect  y=0mm  for  the  images  appearing  in  Figs.  3a  and  3b 
(2mm-diameter  target). 

Fig.4  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and 
without  the  total- variation-m  inimization  (TVM)  enhancement  with  different  contrast 
levels  between  the  target  and  the  background  (case  3).  (a),  without  the  TVM,  1.5:1 
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contrast,  (b),  with  the  TYM,  1.5:1  contrast  .  (c),  absorbed  energy  density  profiles  along 
the  transect  y=0mm  for  the  im  ages  appearing  in  Figs,  la  and  lb  (2:1  contrast).  (d), 
absorbed  energy  density  profile  s  along  the  transect  y=0mm  for  the  images  appearing  in 
Figs.  4a  and  4b  (1.5:1  contrast). 

Fig.5  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and 
without  the  total-variation-  minimization  (TVM)  enhancem  ent  for  three  targets  having 
different  shapes  (case  4).  (a),  exact  image,  (b),  without  the  TVM.  (c),  with  the  TVM  .  (d), 
the  absorbed  energy  density  profiles  along  the  transect  y=0mm. 

Fig.6  Reconstructed  absorbed  energy  de  nsity  images  from  the  three  phantom 
experiments,  (a),  case  1  w  ithout  the  TVM.  (b),  case  2  without  the  TVM,  (c),  case  3 
without  the  TVM.  (d),  case  1  with  the  TVM.  (e),  case  2  with  the  TVM.  (f),  case  3  with 
the  TYM. 

Fig.  7  Recovered  absorbed  energy  density  profiles  along  (a)  y=-7.0mm  crossing  the 
3mm-diameter  target  for  experimental  case  1,  (b)  y=8.0mm  crossing  the  2mm-diameter 
target  for  experimental  case  1,  (c)  y=1.0mm  for  experimental  case  2,  and  (d)  y=6.5mm 
for  experimental  case  3. 

Fig.  8  Recovered  absorbed  energy  density  im  ages  from  the  rat  brain,  (a),  without  the 
TVM.  (b),  with  the  TVM. 
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